O modelo de Solow assume taxa de poupança consante, crescimento da população constante, e progresso tecnológico como exógeno. Um bem homogêneo é consumido e investido em capital $K$ tem técnologia de produção Cobb-Douglas \begin{equation} Y(t) = K(t)^\alpha ( A(t) L(t) )^{1-\alpha}, \quad 0 < \alpha < 1 \tag{1} \end{equation}
Suponha que o estoque de tecnologia $A(t)$ e o estoque de trabalho $L(t)$ ambos cresçam a taxas constantes: \begin{align} &\frac{dA(t)}{dt} = g, \quad &A(t) = A(t)\exp(gt); \ &\frac{dL(t)}{dt} = n, \quad &L(t) = L(t)\exp(nt). \end{align}
Podemos escrever (1) em termos de unidades efetivas de trabalho:
\begin{align} \frac{Y(t)}{A(t)L(t)} &= \frac{1}{(A(t)L(t))^\alpha (A(t)L(t))^{1-\alpha} } K(t)^\alpha (A(t)L(t))^{1-\alpha} \\ y(t) &= k(t)^\alpha \tag{2} \end{align}Letras minúsculas denotam a variável em unidades de trabalho efetivo.
Definiremos a equação de movimento do capital efetivo por trabalhador (que descreve a trajetória da variação do estoque de capital, e, no nosso caso também do produto) como \begin{equation} \frac{dk(t)}{dt} = s y(t) - (n + \delta + g)k(t) \tag{3} \end{equation} \begin{equation} \frac{dk(t)}{dt} = s k(t)^\alpha - (n + \delta + g)k(t) \tag{4} \end{equation}
No estado-estacionário (steady-state) do modelo de Solow temos que $\frac{dk(t)}{dt}=0$. Note que isso implica $\frac{dy(t)}{dt}=0$ no Steady-State.
α = 0.33 # proporção do capital na renda
n = 0.02 # taxa de crescimento da população
δ = 0.01 # taxa de depreciação do capital físico
g = 0.04 # taxa de crescimento do produto
plot(plt_a, plt_b, grid=true, layout=@layout([a;b]))
Igualando o lado direito da equação (4) a zero, temos que o valor para o qual o estoque de capital por unidade efetiva de trabalho converge no estado estacionário:
\begin{equation} k^* = \left( \frac{s}{n+g+\delta} \right)^{1/(1-\alpha)} \tag{5} \end{equation}Tudo mais constante, $k^*$ é maior conforme a taxa de poupança cresce (mas cada vez menos!).
Agora vamos manipular equação (1) de modo a obter uma forma funcional que possa ser estimada por MQO Substiremos $k^*$ na função de produção (1), e dividiremos pelo estoque de trabalho $L(t)$, mas não pelo estoque de tecnologia $A(t)$ (por quê?) \begin{equation} Y(t) = K(t)^\alpha \left( A(t)L(t) \right)^{1-\alpha} \end{equation} \begin{equation} \frac{Y(t)}{L(t)} = \left( \frac{K(t)}{L(t)} \right)^\alpha \left( \frac{A(t)L(t)}{L(t)} \right)^{1-\alpha} \end{equation} \begin{equation} \frac{Y(t)}{L(t)} = \left( A(t)k(t) \right)^\alpha A(t)^{1-\alpha} \end{equation} \begin{equation} \frac{Y(t)}{L(t)} = (k^*)^\alpha A(t)^1 \end{equation} \begin{equation} \frac{Y(t)}{L(t)} = \left( \frac{s}{n+g+\delta} \right)^{\frac{\alpha}{1-\alpha}}A(0)\exp(gt) \end{equation} Tomando logaritmos: \begin{equation} \log \left( \frac{Y(t)}{L(t)} \right) = \log A(0) + gt + \frac{\alpha}{1-\alpha}\ln (s) - \frac{\alpha}{1-\alpha} \ln (n + g + \delta). \end{equation} Para contrastar a equação (6) com os dados temos que tratar $A(0)$. Não temos uma medida concreta da tecnologia inicial em cada país, então seguindo Mankiw et al. [1992] iremos fazer a seguinte suposição: \begin{equation} \log A(0) = \bar{a} + \epsilon_i, \end{equation} em $\bar{a}$ é uma constante e que $\epsilon_i$ é um choque país-específico. Finalmente, temos a equação (7) de Mankiw et al. [1992]: \begin{equation} \log \left( \frac{Y(t)}{L(t)} \right) = \bar{a} + gt + \frac{\alpha}{1-\alpha}\ln (s_i) - \frac{\alpha}{1-\alpha} \ln (n_i + g + \delta) + \epsilon_i. \tag{FE} \end{equation} Suponha $t=0$. Isto é: nosso "um" período $t=0$ vai de 1970 a 2014.
(PA). A renda per capita de estado depende positivamente da taxa de poupança e negativamente da taxa de cresimento populacional. (previsões qualitativas)
(PB). Porque o modelo de Solow assume que os fatores de produção são remunerados de acordo com seus produtos marginais, e só temos dois fatores de produção; então sabemos que a magnitude da elasticidade-poupança da renda per-capita deve ser $\frac{\alpha}{1-\alpha}$, e que a magnitude da elasticidade-"n" da renda per-capita deve ser $-\frac{\alpha}{1-\alpha}$. (lembrando que $\alpha$ é a contribuição do capital para a renda) (previsão quantitativa).
Nosso objetivo é usar técnicas de regressão linear para testar (PA) e (PB).
df_vars = readtable("solow_vars.csv")[:,2:end];
head(df_vars)
| ISO3 | GDPPOP | IOY | GPOP | SCHOOL | dSCHOOL | g | |
|---|---|---|---|---|---|---|---|
| 1 | AGO | 7967.85726432466 | 0.339072699613636 | 0.03082154772623 | 1.20071685034091 | 0.00815095183746249 | 0.0419973074459723 |
| 2 | ALB | 10664.4194020368 | 0.168105318 | 0.00632000422495465 | 2.41369298097727 | 0.0142636035341021 | 0.0329134983094553 |
| 3 | ARE | 64397.5125187937 | 0.364501211477273 | 0.0813291301815251 | 2.15045196909091 | 0.0155644636347753 | 0.0502477150945378 |
| 4 | ARG | 20221.829321369 | 0.16118999175 | 0.0131983404266937 | 2.53402402720455 | 0.00790888325521768 | 0.0503095429744743 |
| 5 | AUS | 43070.6193832596 | 0.280733800045455 | 0.0136215517254629 | 3.35598986259091 | 0.00352776157512386 | 0.0322031572730742 |
| 6 | AUT | 47744.4608470954 | 0.283952697068182 | 0.00282266357823816 | 2.99514657795455 | 0.00542856285494992 | 0.0306597636384636 |
plot(hst1,hst2,hst3,hst4)
plot(lhst1,lhst2,lhst3,lhst4)
sctr1
sctr2
sctr3
df_vars_all = DataFrame()
df_vars_all[:y] = Array(log.(df_vars[:GDPPOP]))
df_vars_all[:ioy] = Array(log.(df_vars[:IOY]))
df_vars_all[:n] = Array(log.(df_vars[:GPOP]+0.05))
df_vars_all[:ioy_minus_n] = df_vars_all[:ioy]-df_vars_all[:n]
df_vars_all[:school] = log.(df_vars[:SCHOOL])
df_vars_all[:school_minus_n] = df_vars_all[:school]-df_vars_all[:n]
df_vars_all[:dschool] = log.(df_vars[:dSCHOOL]);
Forma Funcional 1: \begin{equation} \log \frac{Y_i}{L_i} = \beta_0 + \beta_1 \log s_i + \beta_2 \log(n_i + g + \delta) + \epsilon_i \end{equation}
# Estimação
y = Array(df_vars_all[:y])
X = Array([ones(y) df_vars_all[:ioy] df_vars_all[:n]]) # Primeira Especificação!
N = length(y) # tamanho amostral
p = size(X,2) # número de parâmetros estimados
β = inv(X'X)X'y
3-element Array{Float64,1}:
7.75598
1.98036
-1.79012
ϵ = y - X*β # resíduos
σ2_ϵ = ϵ'ϵ/(N-p) # variância dos resíduos
Σ2_β = σ2_ϵ*inv(X'X) # matriz de covariância do estimador β
3×3 Array{Float64,2}:
2.06777 0.143904 0.678709
0.143904 0.051206 0.0233547
0.678709 0.0233547 0.237586
Uma maneira de testar (PA) é um teste t individual para os parâmetros $\beta_1$ e $\beta_2$
$H_0$: $\beta_i = 0$
$H_1$: $\beta_i \neq 0$
Sabemos que $\frac{\beta_i}{SE(\beta_i)} \sim TDist_{(N-p)}$. Portanto podemos calcular as estatísticas de teste:
t_obs = β./sqrt.(diag(Σ2_β))
3-element Array{Float64,1}:
5.39369
8.75154
-3.67258
pvals = ccdf(TDist(N-p),abs.(t_obs))*2 #*2 porque Normal() é simétrica, e nosso teste é bicaudal.
3-element Array{Float64,1}:
4.7875e-7
6.16518e-14
0.000391292
TSS = (y - mean(y))'*(y - mean(y)) #soma de quadrados totais
RSS = ϵ'ϵ # soma de quadrado dos resíduos
R2 = 1-RSS/TSS # coeficiente de determinação
adj_R2 = 1-(1-R2)*(N-1)/(N-p) # R² ajustado
df_out = DataFrame(); df_out[:coefs] = β; df_out[:t_obs] = t_obs; df_out[:pvals] = pvals; df_out[:R2]= R2;
df_out[:adj_R2] = adj_R2
print(df_out)
3×5 DataFrames.DataFrame │ Row │ coefs │ t_obs │ pvals │ R2 │ adj_R2 │ ├─────┼──────────┼──────────┼─────────────┼──────────┼──────────┤ │ 1 │ 7.75598 │ 5.39369 │ 4.7875e-7 │ 0.525551 │ 0.515868 │ │ 2 │ 1.98036 │ 8.75154 │ 6.16518e-14 │ 0.525551 │ 0.515868 │ │ 3 │ -1.79012 │ -3.67258 │ 0.000391292 │ 0.525551 │ 0.515868 │
reg1 = lm(@formula(y~ioy+n), df_vars_all)
DataFrames.DataFrameRegressionModel{GLM.LinearModel{GLM.LmResp{Array{Float64,1}},GLM.DensePredChol{Float64,Base.LinAlg.Cholesky{Float64,Array{Float64,2}}}},Array{Float64,2}}
Formula: y ~ 1 + ioy + n
Coefficients:
Estimate Std.Error t value Pr(>|t|)
(Intercept) 7.75598 1.43797 5.39369 <1e-6
ioy 1.98036 0.226287 8.75154 <1e-13
n -1.79012 0.487428 -3.67258 0.0004
Forma Funcional 2: \begin{equation} \log \frac{Y_i}{L_i} = \gamma_0 + \gamma_1 [\log s_i - \log(n_i + g + \delta)] + \eta_i \end{equation}
reg2 = lm(@formula(y~ioy_minus_n), df_vars_all)
DataFrames.DataFrameRegressionModel{GLM.LinearModel{GLM.LmResp{Array{Float64,1}},GLM.DensePredChol{Float64,Base.LinAlg.Cholesky{Float64,Array{Float64,2}}}},Array{Float64,2}}
Formula: y ~ 1 + ioy_minus_n
Coefficients:
Estimate Std.Error t value Pr(>|t|)
(Intercept) 7.28953 0.224388 32.4863 <1e-53
ioy_minus_n 1.93808 0.185267 10.461 <1e-16
# o "implied α" de Mankiw et al. [1992] na Tabela 1 vem da Forma Funcional 2 e da Forma Estrutural (FE).
β₁ = coef(reg2)[2]
implied_α = β₁/(1+β₁)
Vamos testar agora a previsão (PB): que os coeficientes de $log (s)$ e $log(n +g+\delta)$ têm mesma magnitude, mas sinais contrários. Para tal compararemos a acurácia dos dois modelos que já estimamos. Primeiro o modelo irrestrito (UR), no qual o parâmetros que multiplicam o log da taxa de poupança $log(s)$ e $log(n+g+\delta)$ podem ser diferentes. (nossa Forma Funcional 1!) \begin{equation} \log \frac{Y_i}{L_i} = \beta_0 + \beta_1 \log s_i + \beta_2 \log(n_i + g + \delta) + \epsilon_i \end{equation} Então o modelo restrito (R), no qual fixamos que os parâmetros têm mesma magnitude mas sinais contrários. (nossa Forma Funcional 2!) \begin{equation} \log \frac{Y_i}{L_i} = \gamma_0 + \gamma_1 [\log s_i - \log(n_i + g + \delta)] + \eta_i \end{equation}
$H_0$: $\beta_1$ = $\beta_2$ = $\gamma_1$.
$H_1$: caso contrário.
\begin{equation} F_{obs} = \frac{SSR_R - SSR_{UR}}{SSR_{UR}} \frac{N - q_{UR}}{q_R - q_{UR}} \sim F(q_{UR}-q_{R}, N-q_{UR}) \end{equation}No nosso caso: $q_{UR} = 3$, $q_{R} = 2$, $N = 101$.
reg1 = lm(@formula(y ~ ioy + n), df_vars_all) # estimamos o modelo conforme a Forma Funcional 1
reg2 = lm(@formula(y ~ ioy_minus_n), df_vars_all) # estimamos o modelo conforme a Forma Funcional 2
resid1 = y - predict(reg1) # calculamos os resíduos da FF1
resid2 = y - predict(reg2) # calcularmos os resíduos da FF2
SSR_R = resid2'resid2 # Soma do quadrado dos resíduos da FF1
SSR_UR = resid1'resid1 # Soma do quadrado dos resíduos da FF2
q_UR = length(coef(reg1)) # Número de parâmetros na FF1
q_R = length(coef(reg2)) # Número de parâmetros na FF2
N = length(y) # Tamanho amostral
F_obs = ((SSR_R-SSR_UR)/(q_R-q_UR))/(SSR_UR/(N-q_UR)) # Estatística do teste F que elaboramos
pval_F = ccdf(FDist(q_UR-q_R, N-q_UR), F_obs)
# pval = 1 indicando que a restrição de que β₁ = β₂ não piora significativamente o "fit" do modelo.